On Quantum Effects near the Liquid- Vapor Transition in Helium 
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The liquid-vapor transition in 3 He and 4 He is investigated by means of path-integral molecular 
dynamics and the quantum virial expansion. Both methods are applied to the critical isobar and 
the critical isochore. While previous path-integral simulations have mainly considered the lambda 
transition and superfluid regime in 4 He, we focus on the vicinity of the critical point and obtain 
good agreement with experimental results for the molar volume and the internal energy down to 
subcritical temperatures. We find that an effective classical potential that properly describes the 
two-particle radial distribution function exhibits a strong temperature dependence near the critical 
temperature. This contrasts with the behavior of essentially classical systems like xenon, where the 
effective potential is independent of temperature. It is conjectured that, owing to this difference in 
behavior between classical and quantum-mechanical systems, the crossover behavior observed for 
helium in the vicinity of the critical point differs qualitatively from that of other simple liquids. 



I. INTRODUCTION 

Crossover phenomena have enjoyed a renewed atten- 
tion in recent years, both from the theoretical (see, e.g., 
Refs. and from the experimental side (cf. Refs. [(] |g); 
see also Ref. ^ and references therein. This concerns 
in particular the crossover from mean-field-like to Ising- 
like critical behavior upon approach of the critical point. 
The accurate numerical determination of crossover scal- 
ing functions for the isothermal compressibility and the 
liquid-vapor coexistence curveO has motivated-the reex- 
amination of experimental data for 3 He and XeJliJ Xenon, 
with a very high molar mass and a relatively high critical 
temperature, is essentially a classical system, while the 
critical point for 3 He occurs at a temperature and density 
where quantum effects are expected to be non- negligible. 
This is also expressed by the de Boer parameter A* $3 For 
monatomic gases of atomic mass m, that are described 
by a Lennard-Jones potential with parameters s and a, 
this parameter is defined as A* = h/ ay/me. For 3 He its. 
value is A* = 3.08, compared to A* = 0.064 for xenon.t 3 ] 
Nevertheless, the nature of the critical point itself is the 
same for both fluids, because the critical fluctuations 
dominate over the quantum-mechanical fluctuations at 
temperatures sufficiently close to the critical tempera- 
ture T c . Thus, the values of the critical exponents are 
not affected. For the crossover region, the situation is 
less clear-cut. Since the correlation length now has a fi- 
nite (although large) value, nonuniversal behavior may 
be expected for different systems and has actually been 
observed.Q However, according to theory, this nonuniver- 
sality is largely determined by the so-called cut-off pa- 
rameter (the wavenumber corresponding to some micro- 
scopic characteristic length) which in most simple fluids 
has a very similar valueJH Furthermore, also in numerical 
simulations of Ising-like systems a high degree of univefc 
sality has been observed for crossover scaling functions.tll 
Thus, it came as quite a surprise that the crossover be- 
havior for 3 He exhibits a marked difference from that of 



xenon. Owing to the short range of the interactions, 
these systems cannot be expected to complete the full 
crossover to mean-field-like behavior before leaving the 
critical region [where t <C 1, with t = (T—T c )/T c ]. Never- 
theless, for T > T c , the crossover behavior of the isother- 
mal compressibility of xenon turned out to agree very well 
with the numerical data for the three=dimensional Ising 
model with varying interaction rangeji9 whereas the cor- 
responding experimental data for 3 He seemed to be de- 
scribed by a qualitatively different curve. For t > 0.01, 
the compressibility appeared essentially suppressed com- 
pared to the crossover scaling function. In Ref. 11, this 
difference in behavior was conjectured to be related to 
quantum-mechanical effects: The critical compressibility 
of 3 He would be enhanced due to quantum fluctuations, 
which are temperature dependent. Hence, this contribu- 
tion is expected to decrease appreciably within the (high- 
temperature) crossover region, effectively leading to an 
additional reduction of the compressibility upon increase 
of the temperature. Clearly, it is only the compressibility 
due to the thermal fluctuations which is described by the 
various theoretical expressions for the crossover scaling 
function. 

In order to gauge the quantum-mechanical contribu- 
tion to the compressibility, a comprehensive theoretical 
description of the critical behavior, including the role of 
quantum effects, is required; such a description should 
encompass the temperature dependence of this contribu- 
tion in the vicinity of T c . Indeed, for a weakly inter- 
acting Bose fluid a scaling function has been calculated 
describing the crossover from criticalLty (i.e., the lambda 
point) to ideal Bose-gas behaviorMEil In particular, it 
was found that a mapping of the Hamiltonian for the 
Bose gas onto a classical spin model yields a Landau- 
Ginzburg- Wilson Hamiltonian with a quartic tepu that 
has a strongly temperature-dependent coefficient .E£l Since 
this coefficient plays a pivotal role in crossover scaling 
functions, a corresponding effect in the vicinity of the 
liquid-vapor critical point would definitely affect the na- 
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ture of the crossover from Ising-like to mean-field-like 
critical behavior. These considerations have motivated 
us to examine the behavior of helium close to the critical 
point by means of quantum-mechanical numerical meth- 
ods. In particular, it is of interest to see whether the 
magnitude of quantum effects indeed changes apprecia- 
bly over the crossover region, as conjectured in Ref. [fl|. 
In the present study, we pay some attention to 3 He, but 
our main focus is on 4 He. Due to the higher mass of 
He, quantum effects are less pronounced, which is only 
reinforced by the correspondingly higher critical temper- 
ature. This facilitates the numerical calculations con- 
siderably, as will be outlined in Sec. [I A. At the same 



time, the quantum effects are still clearly visible in the 
crossover region and the isothermal compressibility in- 
deed exhibits a deviation from the predicted crossover 
curve, similar to that found for 3 He.Ej For comparison, 
4 He has a de Boer parameter A* = 2.67, compared to the 
above-mentioned value A* = 3.08 for 3 He.E3 

For completeness we remark that 4 He has been ex- 
tensively studied by means of the path-integral Monte 
Carlo method, especially in the context of the lambda 
transition taking place at 2.17 K, cf. Ref. [l7[ Also the 
high-density region has been explored in this wayjla but 
we are not aware of numerical studies in the vicinity of 
the liquid- vapor critical point. 

The outline of this paper is as follows. In Sec. || wc in- 
troduce the methods that we have appli ed: p ath-integral 
molecular dynamics is discussed in Sec. [I A , followed by 
the quantum virial expansion in Sec. [IB. Section III 



contains all our main results, namely data for the atomic 
volume of 4 He and for its kinetic, potential, and inter- 
nal energy along the critical isobar, as well as for the 
kinetic and potential energy of 4 He and 3 He along their 
respective critical isochores. Furthermore, results for the 
effective pair potential of 3 He in the vicinity of the critical 
temperature are presented. Our conclusions are summa- 
rized in Sec. 



II. METHODS 

A. Path-integral Molecular Dynamics 

Path-integral Monte Carlo (PIMC)0S and path- 
integral molecular dynamics (PIMD)E!] are well- 
established methods to calculate thermodynamic prop- 
erties of many-particle quantum systems. Path inte- 
gral techniques exploit the possibility to represent the 
quantum-mechanical partition function of point partL 
cles as a classical partition function of closed polymers.EZl 
Neighboring beads in the polymer are coupled via elas- 
tic springs with stiffness k — mM 2 / j3 2 h 2 , where m is 
the mass of the point particle, the so-called Trotter num- 
ber M represents the number of beads in the polymer, 
and (3 = l/k^T. Interaction between polymers only 
takes place between those monomers that have the same 



index within the respective polymers they belong to, 
and the effective temperature in the isomorphic classi- 
cal representation is given by TM . For more details, 
wc rcfci: the , .reader to the original literature and review 
articlesHHBB 

Depending on the nature of the problem under study, 
Monte Carlo methods are preferable to molecular dy- 
namics or vice versa. PIMC has so far been the method 
of choice for finite-temperature simulations of condensed 
helium, in particular when exchange effects played an 
important role. Far away from the superfluid regime, 
however, exchange effects can be neglected, 113 so that 
PIMC is not necessarily advantageous to PIMD for the 
study of the liquid-vapor transition. For example, in the 
gas phase, ballistic trajectories can be realized in PIMD, 
which allows rapid changes in the configuration. In con- 
ventional Monte Carlo methods, motions are constrained 
to be diffusive and large correlation times may be ex- 
pected. Also, in isobaric simulations, volume moves are 
done at no extra cost in a molecular dynamics simula- 
tion, while Monte Carlo methods require the evaluation 
of the total energy of the system. 

Using an appropriate representation of the internal co- 
ordinates of the chaip molecules, it is possible to avoid 
inefficient samplingril; e.g., in the regular representation, 
where each bead in a chain has the same dynamical in- 
ertia, the time-step discretization has to be chosen pro- 
portional to M~ 2 . In this study, we have used a repre- 
sentation of the chain molecules in terms of the center- 
of-mass coordinate and the eigencoordinates of the free 
particles, defined in Eq. (|[) below, which makes it pos- 
sible to work with time-steps independent of M. While 
the "dynamical" center-of-mass coordinate is chosen to 
be identical with the real mass to, different masses m q 
are attributed to each eigenmode q. An efficient choice is 
m q = (k + k q )m/k, where k is an adjustable parameter 
and k q is the stiffness associated with the eigenmode q. k 
is conveniently chosen such that in the condensed phase 
all modes move on approximately the same time scale. 

The main disadvantage of PIMD is that finite-time step 
errors accumulate for certain quantities when the Trotter 
number M becomes large. While the average potential 
energy and the virial estimator for the kinetic energycfJ 
do not suffer from this effect, the so-called primitive es- 
timator for the kinetic energy, i^ pr i m , does. An estima- 
tor is a function whose average value corresponds to the 
expectation value of a property of interest. For an in- 
dividual configuration, however, the association of the 
actual value of the estimator with the actual, value of the 
property is meaningless. K pr i m is given bye 4 ] 

Q h N M 

K plim = -Nk B TM 2 - 2 E E ( r - - r -+i) 2 • W 

i=l t=1 

where iv represents the position of the rth monomer in 
polymer i, and k is the stiffness introduced at the begin- 
ning of this section. In addition to the above-mentioned 
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systematic error, K pr - lul has been shown to have large 
statistical errors fcpJarge Trotter numbers M, even in 
PIMC simulations.cJ However, a simple trick remedies 
both of these shortcomings .E3 If ^Nk^TM 2 is replaced 
by the actual "dynamic" kinetic energy, a new estimator 
can be defined, 



^ N M 

-^prim = - ^ ^ ( m g v ig — k q vf qj 



(2) 



1=1 q=l 



where v' iq denotes the velocity associated with the eigen- 
coordinate 
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2-Kiqr/M 



(3) 



It turns out that the new estimator ifp r i m is superior to 
the virial estimator. In the solid phase, e.g., the virial 
estimator misses the center-of-mass motion of the sim- 
ulation box, leading to a systematic difference between 
the primitive and the virial estimator of ZkgT /2 for the 
whole system. While this effect is negligible for large 
particle numbers, the same shortcoming has more serious 
implications in the gas phase: The virial estimator sys- 
tematically underestimates the kinetic energies because it 
does not include the kinetic energy of the center-of-mass 
motion of independent clusters. In the condensed phase, 
the statistical errors of both estimators are similar .E3 Un- 
fortunately, -Kprim cannot be used in Monte Carlo (MC) 
simulations, because MC is not based on deterministic 
dynamics but solely on stochastic dynamics. It is worth 
mentioning that the average of both estimators tends to 
the exact thermal expectation value of the kinetic energy 
with increasing Trotter number M. For methods to per- 
form PIMD simulations at constant pressure, we refer to 
Refs. HfH 

All simulations, as well as the vifial expansions, are 
based on the Aziz HFD-B potential,E3 which is considered, 
one of the best known interatomic model potentials.tll 
It consists of a Hartree-Fock (exponential) short-range 
repulsive term and algebraic long-range attractive terms 
(1/r 6 , 1/r 8 , 1/r 10 ). The cut-off radius that we used in 
our simulations was r c — 10 A. The particle number in 
all simulations was N = 500 and the Trotter number 
varied between M = 1 for classical simulations and M = 
64 for the quantum-mechanical 4 He simulations at the 
lowest temperatures, keeping TM > 200 K. The usual 
corrections of the order of 1/M 2 were applied to the final 
data.E3 For the 3 He simulations, TM > 350 K was used 
with a maximum value of M = 128. We have used a 
cubic simulation box with periodic boundary conditions; 
the length of each simulation amounted to at least 50 000 
time steps. 

Note that our simulations mostly address the identifi- 
cation of quantum effects in the first derivatives of the 
thermodynamic potential (internal energy, volume, etc.), 
while an estimation of second derivatives such as the com- 
pressibility (discussed in Ref. Ill]) is not attempted: This 



would require much larger system sizes and is not feasible 
with our computer resources. 



B. Virial Expansion 

An alternative method for calculating thermodynamic 
properties of quantum gases is by means of the quantum 
virial expansion. In the original approach (see Ref. |l3| , 
Chapter 6 and references therein), a series expansion in 
h is obtained for every vi rial coefficient. Here, a more 
efficient technique is used,BaEJ which is briefly outlined 
below. 

The first correction to the internal energy of an 
ideal gas arises from the pair interaction U12 (/?, V) of 
(quantum-mechanical) particles. A pair of particles is 
described by a center of mass mode, which can be treated 
classically in the second-order virial expansion, and the 
relative coordinate r, which is confined to a (spherical) 
volume V. Throughout the derivation, finite-volume cor- 
rections jare ignored. For N 2 /2 particle pairs, U(N, V, (3) 
becomes!!] 

37V N 2 ( /iV x 

U(N, V, (3) = — + — u 12 ({3, V)+0lN^ v 



(4) 

where U12 (/?, V) is calculated according to 

un(P,V) = -^5 log / d 3 r (r | e -0&-+«») | r ) 
op Jy 

+^l«g/ d 3 r{r\e-^ | r) . (5) 

where i rc i denotes the operator for the kinetic energy 
associated with the relative motion of two particles and 
v 12 is their potential energy. The two integrands on the 
right-hand side of Eq. @ are the diagonal elements of 
the density matrix p(r,r',/3), 



p(r,r',P) = (r | exp[-/3(t re l + ««)] | r') 



(G) 



and of its noninteracting counterpart. 

Up to a normalization factor, which is irrelevant for 
the calculation of u±2(P, V) in Eq. (||), the radial distri- 
bution function g(r) is given by the diagonal elements of 
p(r, r',/3). Taking into account that J v d 3 rg(r) — > V in 
the thermodynamic limit, it is possible to rewrite Eq. (j^) 
in the more familiar form 



2 dB 2 (f3) 
V d(3 ' 



(7) 



where, like for classical systems, the second virial coeffi- 
cient £?203) can be expressed in terms of g(r): 



B 2 (/3)=-27r/ drr'[g 12 (r)-g (r) 
Jo 



(8) 
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Here 1712 (r) and go (r) denote the radial distribution func- 
tion in the interacting and non-interacting case, respec- 
tively. 

The diagonal elements p(r, r, 0) can be calculated by 
exploiting the semi-group property of the density opera- 
tor 

p(v, r', 2/3) = J d 3 r"p(r, r", 0)p(r", r', (3) . (9) 

Thus we can obtain the low-temperature density matrix 
at temperature T by squaring the density matrix at tem- 
perature 2T. For n iterations, the starting temperature 
has to be chosen as 2™T. For the highest temperature 
of the iteration process, it is possible to use the so-called 
primitive decomposition for p(r, r', /3/M), which under- 
lies the path-integral simulations presented in this pa- 
per as well as most other path- integral simulations. One 
of the advantages of the squaring procedure over path- 
integral simulations is that the required numerical ef- 
fort only scales logarithmically with inverse temperature. 
Hence, it is easy to minimize discretization errors (for 
two-particle systems) using squaring techniques at low 
temperatures. At a given (high) temperature, the sys- 
tematic error is proportional to 1/M 2 , just like in path- 
integral simulations .E3 

For spherically symmetrical potentials, Eq. (^) can be 
reduced to a sum of one-dimensional integrations by de- 
composing the density matrix into |Cantributions belong- 
ing to different angular momentaEauj In practice, the 
squaring is done in terms of simple matrix multiplication 
by discretizing the variable r. Of course, a cut-off r c has 
to be introduced at a reasonably large value of r. This 
induces artificial behavior at the boundary not found in 
an infinitely large system, namely that go (r) tends to 
zero as r approaches r c . Therefore, the integration in 
Eq. (|J) has to be confined to the region where boundary 
effects are negligible. Alternatively, one may normalize 
the integrand in Eq. (j^) by go( r )j which results in a fast 
convergence of B% with r c . 

Quantum effects in the calculation of i?2(/3) will be- 
come important when the thermal wavelength A(/3) — 
hj \f2wmk^T of the free particle is in the order of or 
larger than the distance at which the interatomic poten- 
tial is minimum. We want to illustrate this in the case 
of 4 He at a temperature T = 10 K, where A ~ 2.8 A. In 
Fig. [I], the two-particle radial distribution function g(r) 
is shown for a pair of "classical" helium atoms and a pair 
of 4 He atoms. The maximum in g(r) for the quantum- 
mechanical calculation is shifted by about 0.59 A with re- 
spect to the classical equilibrium distance and the height 
of the maximum (relative to go = 1) is decreased by a 
factor of about 5.8. Thus, an effective classical poten- 
tial Ves(r) that would result in a similar g(r) of 4 He at 
T = 10 K as the quantum-mechanical calculation would 
have a strongly reduced binding energy with respect to 
the original Aziz potential and an equilibrium distance 
shifted by 0.59 A. Note that also the curvature of the 



effective potential would be different from the original 
potential. 

In Fig. |], PIMD simulation results, taken at the 
same temperature and the critical pressure P c ^ — 
0.22746 MPa[], are included as well. The agreement be- 
tween the virial expansion and the simulation is very 
good. For the classical system, small differences in g(r) 
can be seen between both methods that can be attributed 
to three-body effects, which are neglected in the second- 
order virial expansion. 

III. RESULTS 
A. Critical isobar of 4 He 

As a first test of the PIMD simulation and the virial 
expansion, we have calculated the atomic volume as a 
function of temperature on the critical isobar, P Cj 4 = 
0.22746 MPa.O While for 3 He the diameter-^ = (pn q + 
Pvapor)/2 has a slope that is almost zero,E3 for 4 He it 
has a clearly positive slope. Thus, while the isobaric 
thermal expansion coefficient V~ 1 (dV/dT)p only has a 
finite peak as a function of temperature for P > P c , 
it diverges at the liquid-vapor transition temperature 
for P < P c . The results for classical and quantum- 
mechanical calculations are presented in Fig. ^ and com- 
parison is made to a phenomenological widft-jjange equa- 
tion of state based on experimental data.EilH As can be 
seen, the overall agreement between the experimental re- 
sults and the quantum-mechanical simulations is very 
good. This demonstrates both the quality of the Aziz 
potential and the fact that the (neglected) three-body 
effects play a remarkably small role here. Only near the 
critical temperature, T C) 4 = 5.1953 K,lJ small systematic 
deviations occur, which however may be tentatively at- 
tributed to finite-size effects in the simulations. In addi- 
tion, the performance of the phenomenological equation 
of state itself might deteriorate in the critical region. The 
second-order quantum virial expansion, included in the 
same figure, exhibits deviations from the experimental 
curve below a temperature of approximately 7 K. 

For comparison, data for "classical helium" are in- 
cluded in Fig. ^ as well. The absence of quantum fluctu- 
ations leads one to expect an increase in both the crit- 
ical temperature and the critical pressure. Thus, the 
classical system should undergo a first-order phase tran- 
sition at the critical pressure P c ,i, as is indeed borne 
out by the numerical data. The corresponding transi- 
tion temperature T\ could be located at approximately 
10 K; here both the fluid and the gas phase were sta- 
ble for the duration of the simulation. Hysteresis effects 



"The subscript '4' will generally be used to indicate quanti- 
ties pertaining to 4 He, whereas a subscript '3' refers to 3 He. 
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and, possibly, slowing down owing to the vicinity of the 
critical point, may have affected the accuracy of this es- 
timate. Since neither the pressure P Cj 4 nor the corre- 
sponding transition temperature play a particular role 
for the classical system, no further attempts were made 
to improve the estimate of T\. We note only that the 
critical temperature X^ lass > T c ^ of the classical sys- 
tem will be even higher than T\ . The magnitude of the 
shift of the transition temperature is a clear indication 
of the importance of quantum-mechanical effects in the 
vicinity of the critical pejnt. Indeed, from the principle 
of corresponding statestia one finds (approximating the 
interaction potential by a Lcnnard- Jones potential with 
a = 2.56 A and e/k B = 10.22 K, Ref. \H§ T c class = 13 K 
and F c class = 1.1 MPa. 

An interesting effect can be observed in the kinetic en- 
ergy Tki rl for 4 He, shown in Fig. |^. Upon lowering the 
temperature, Tkj n suddenly rises in the vicinity of the 
critical temperature. This rise is purely related to the 
increase in the density of the system. It should be noted 
that the behavior of Tki n close to T c is likely to be affected 
by finite-size effects, as these effects generally shift the 
top of the coexistence curve toward temperatures above 
the true critical temperature. The break-down of the 
second-order virial expansion due to this increase in den- 
sity occurs already slightly above T Cj 4, namely around 
5.5 K. The potential energy V po t, depicted in the same fig- 
ure, increases with temperature, as usual. We note, how- 
ever, that this tendency of decreasing kinetic energy and 
increasing potential energy upon increase of the temper- 
ature near T c is not a necessity: A quantum-mechanical 
model for molecular ordering of rotors on a surface rather 
showed the opposite trend, where (Tki n ) increased and. 
(Vpot) decreased near T c upon increasing temperature.Lj 

Figure ^ shows the internal energy alp«|g-|the critical 
isobar, both as obtained from experiment! 31 !! 33 ! and as cal- 
culated by means of PIMD. Although the experimen- 
tal data lie systematically above the numerical ones for 
T > 6 K, the overall agreement is certainly appreciable. 

Finally, we mention that the average classical potential 
energy (Upot), which is not shown in Fig. |[ exhibits a 
clear jump, as expected for a first-order transition. 



B. Critical isochore of He 

We now-Aum to the critical isochore, p Cj 4 = 0.017399 
mol/cm 3 .!?_!] Figure |5| shows both the kinetic and the po- 
tential energy per atom, along with the results of the 
quantum virial expansion. As can be seen, the agreement 
is remarkably good for the kinetic energy, even in the crit- 
ical region. For the potential energy, on the other hand, 
the agreement is not so good, even at relatively high tem- 
peratures. In the same figure, we have also included the 
potential energy for "classical helium" as obtained from 
a second-order and the third-order virial expansion: Here 
the agreement is much better. However, a remark on the 



density is in order here. As mentioned before, quantum 
fluctuations generally lead to a lower critical tempera- 
ture and hence, at the same pressure, to a higher den- 
sity. However, the critical pressure decreases, and the net 
effect is that fluids in which quantum-mechanical effects 
are nonnegligible have a lower critical density than would, 
be expected from the principle of corresponding states.EJ 
This is also nicely illustrated by the critical properties 
of 3 He, which is basically described by the same pair 
potential as 4 He, but has an even lower critical tempera- 
ture, pressure, and density, entirely due to its lpwer mass 
and consequentially larger de Boer parameter .t2l Indeed, 
a quantum-mechanical version of the principle of corre- 
sponding states can be formulated, in which all devia- 
tions froni|-Glassical behavior are parametrized by this 
parameter.!! 3 ! From this slight digression, we conclude 
that the classical data, as shown in the figure, pertain 
to an isochore that is (for the classical system) a subcrit- 
ical one, which might (in addition to the higher order of 
the virial expansion) explain the quite reasonable agree- 
ment. One remaining point, then, is that this isochore 
must cross the vapor branch of the coexistence curve at 
an unknown temperature. 

Just like along the critical isobar, the PIMD results 
for the internal energy along the critical isochore (Fig. ^) 
exhibit good agreement with the experimental data. 



C. Critical isochore of ' He 

Finally, we consider 3 He on its critical isochore, p c ^ = 
0.01374 mol/cm 3 .a Figure ^| shows both the kinetic and 
the potential energy per atom as obtained from PIMD 
simulations and from the second-order quantum virial ex- 
pansion. As for 4 He, the agreement between both types 
of calculations is very good for the kinetic energy and 
rather poor for the potential energy. Also the overall be- 
havior of both energies is similar to that found for He, 
except that in Fig. |] one cannot observe the formation of 
a "plateau" in (Tki n ) at lower temperatures. This is pre- 
sumably due to the fact that the simulations for 3 He do 
not quite reach the critical temperature, T c ^ — 3.317 K 
(see Ref. [32| the value was converted to the Tgo tem- 
perature scale here) , where the flattening of the curve is 
expected to set in. 

In principle, one might use the kinetic energy to de- 
fine an effective, nonlinear temperature scale in which 
the role of the quantum fluctuations has been taken into 
account. By studying the crossover scaling function for, 
e.g., the compressibility on such a redefined temperature 
scale, one could examine the role of quantum effects in 
the deviations observed in Ref. [ll]. However, the defini- 
tion of such a temperature scale requires a very accurate 
knowledge of the kinetic energy, in particular for temper- 
atures very close to T c , since crossover scaling functions 
are studied on a logarithmic scale in the reduced tem- 
perature t. Unfortunately, the numerical accuracy of our 
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PIMD data did not warrant a meaningful, direct reex- 
amination of the crossover scaling functions. As an al- 
ternative, we provide here qualitative evidence justifying 
the conjecture of Ref. |ll| that the influence of quantum 
effects changes appreciably within the crossover regime. 
To this end, we have considered the effective potential 
V e s as defined by 

where g(r) is the two-particle correlation function. The 
latter quantity, in turn, can be obtained from the 
virial expansion. Both V e s and g(r) are temperature- 
dependent quantities. The resulting effective potential is 
shown in Fig. for three different temperatures, namely 
T = 3.125 K (just below T c ), T = 6.125 K (roughly twice 
T c ), and T = 12.125 K (somewhat above the highest 
temperature studied in Ref. |Il]). One observes that the 
effective potential changes rather dramatically under this 
temperature variation: The depth of the potential well 
decreases by roughly a factor of 3 when the temperature 
is reduced from 12.125 to 3.125 K, and its width increases 
accordingly. This is in concordance with the observed en- 
hancement of the compressibility over this temperature 
region. 

IV. DISCUSSION AND CONCLUSION 

In this paper, we have presented the results of a nu- 
merical study of 4 He and 3 He in the vicinity of their re- 
spective liquid-vapor critical points. Both path-integral 
molecular dynamics and quantum virial expansions have 
been applied. I— . 

For 4 He, we have shown that the Aziz potentiaEJ em- 
ployed in the PIMD calculations-yields very good agree- 
ment with experimental resultsEU'Eij for the atomic vol- 
ume at the critical isobar over a wide temperature range, 
including temperatures that lie considerably below the 
critical temperature. Equally good agreement is found 
for the internal energy. The second-order quantum virial 
expansion shows good agreement for the atomic volume 
down to roughly two Kelvin above the critical tempera- 
ture. The simulations also demonstrate the importance 
of quantum effects in 4 He in the vicinity of its critical 
point: A classical system with the same potential has a 
critical temperature that is more than twice as high. 

In addition, we have studied the kinetic and potential 
energy of 4 He and 3 He along their critical isochores. For 
the kinetic energy, there is good agreement between the 
quantum virial expansion and the simulation results in 
both cases, down to quite low temperatures, but for the 
potential energy the agreement is not so good. Higher- 
order terms in the virial expansion would probably allevi- 
ate the discrepancy, but are tedious to calculate. On the 
other hand, simulation results for the isochoric internal 
energy of 4 He agree well with experimental results. 



One of our original goals was the definition of an effec- 
tive, nonlinear temperature scale based on the kinetic 
energy. Such a temperature scale might prove useful 
for the reexamination of the crossover scaling functions 
of helium, which were experimentally observed to differ 
qualitatively from their counterparts for other simple liq- 
uids; these deviations have been conjectured to originate 
from quantum-mechanical effects]! 1 ] Indeed, since (part 
of) these effects would now be incorporated in the new 
temperature scale, the crossover scaling functions should 
exhibit a better agreement with the (classical) predic- 
tion if the conjecture were correct. However, it is very 
difficult to reach the required numerical accuracy in the 
kinetic energy; in addition, very large particle numbers 
are needed in the immediate vicinity of the critical point, 
in order to circumvent finite-size effects. Thus, we have 
opted for a different strategy, namely the calculation of 
the effective pair potential for 3 He by means of the quan- 
tum virial expansion evaluated on the critical isochorc. 
We find that this pair potential exhibits a strong varia- 
tion with temperature in the crossover region. This is in 
good accord with the conjecture of Ref. O: Namely, not 
only are there strong quantum effects in^He close to its 
liquid-vapor critical point, as was already evident from 
the values of the critical amplitudes (cf., e.g., Ref. [32]), 
but these quantum effects also vary considerably within 
the crossover region. Such a variation might then ex- 
plain the observed apparent depression (relative to the 
predicted crossover scaling function) of the compressibil- 
ity above T c when the temperature is increased. 

ACKNOWLEDGMENTS 

E.L. acknowledges stimulating comments by Michael 
E. Fisher, as well as support from the National Science 
Foundation (through Grant No. CHE 99-81772 to M. E. 
Fisher). M.H.M. thanks Kurt Binder for useful discus- 
sions and acknowledges support from the BMBF through 
Grant 03N6015 and from the Materialwissenschaftliche 
Forschungszentrum. 



1 C. Bagnuls and C. Bervillier, Phys. Rev. B 32, 7209 (1985). 

2 M. A. Anisimov, S. B. Kiselev, J. V. Sengers, and S. Tang, 
Physica A 188, 487 (1992). 

3 E. Luijten, H. W. J. Blote, and K. Binder, Phys. Rev. Lett. 
79, 561 (1997); Phys. Rev. E 56, 6540 (1997). 

4 S. Caracciolo, M. S. Causo, A. Pelissetto, P. Rossi, and E. 
Vicari, Nucl. Phys. Proc. Suppl. 73, 757 (1999). 

5 E. Luijten, Phys. Rev. E 59, 4997 (1999). 

6 G. Meier, D. Schwann, K. Mortensen, and S. Janfien, Eu- 
rophys. Lett. 22, 577 (1993). 



G 



7 M. A. Anisimov, A. A. Povodyrev, V. D. Kulikov, and J. 
V. Sengers, Phys. Rev. Lett. 75, 3146 (1995). 

8 Y. B. Melnichenko, M. A. Anisimov, A. A. Povodyrev, 
G. D. Wignall, J. V. Sengers, and W. A. Van Hook, Phys. 
Rev. Lett. 79, 5266 (1997). 

9 M. A. Anisimov and J. V. Sengers, in Equations of State 
for Fluids and Fluids Mixtures, edited by J. V. Sengers, 
R. F. Kayser, C. J. Peters, and H. J. White, Jr. (Elsevier, 
Amsterdam, 2000). 

10 E. Luijten and K. Binder, Phys. Rev. E 58, R4060 (1998); 
59, 7254(E) (1999). 

11 E. Luijten and H. Meyer, Phys. Rev. E 62, 3257 (2000). 

12 J. de Boer, Physica 14, 139 (1948). 

13 J. O. Hirschfelder, C. F. Curtiss, and R. Byron Bird, Molec- 
ular Theory of Gases and Liquids (Wiley, New York, 1954). 

14 E. Luijten and K. Binder, Europhys. Lett. 47, 311 (1999). 

15 M. Rasolt, M. J. Stephen, M. E. Fisher, and P. B. Weich- 
man, Phys. Rev. Lett. 53, 798 (1984). 

16 P. B. Weichman, M. Rasolt, M. E. Fisher, and M. J. 
Stephen, Phys. Rev. B 33, 4632 (1986). 

17 D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995). 

18 D. M. Ceperley, R. O. Simmons, and R. C. Blasdell, Phys. 
Rev. Lett. 77, 115 (1996). 

19 D. M. Ceperley and M. H. Kalos, in Monte Carlo Methods 
in Statistical Physics, Topics Curr. Phys. Vol. 7, edited by 
K. Binder (Springer, Berlin, 1979). 

20 K. E. Schmidt and D. M. Ceperley, in The Monte Carlo 
Method in Condensed Matter Physics, Topics Appl. Phys. 
Vol. 71, edited by K. Binder (Springer, Berlin, 1992). 

21 M. E. Tuckerman, B. J. Berne, G. J. Martyna, and M. L. 
Klein, J. Chem. Phys. 99, 2796 (1993). 

22 R. P. Feynman and A. R. Hibbs, Quantum Mechanics and 
Path Integrals (McGraw-Hill, New York, 1965). 

23 D. Marx and M. H. Miiser, J. Phys.: Condens. Matter 11, 
R117 (1999). 

24 M. Herman, E. F. Bruskin, and B. J. Berne, J. Chem. Phys. 
76, 5150 (1982). 



25 P. Schoffel and M. H. Miiser, |cond-mat/0009379 . 

26 G. J. Martyna, A. Hughes, and M. E. Tuckerman, J. Chem. 
Phys. 110, 3275 (1999). 

27 R. A. Aziz, F. R. W. Court, and C. C. K. Wong, Mol. Phys. 
61, 1487 (1987). 

28 M. Suzuki, in Quantum Monte Carlo Methods in Equilib- 
rium and Nonequilibrium Systems, edited by M. Suzuki 
(Springer, Berlin, 1987). 

29 R. G. Storer, J. Math. Phys. 9, 964 (1968). 

30 D. Thirumalai, E. J. Bruskin, and B. J. Berne, J. Chem. 
Phys. 79, 5063 (1983). 

31 E. W. Lemmon, M. O. McLinden, and D. G. Friend, 
in NIST Chemistry WebBook, NIST Standard Reference 
Database Number 69, edited by W. G. Mallard and P. J. 
Linstrom (National Instit ute of Standards and Tec hnology, 
Gaithersburg, 2000). See http://webbook.nist.gov. 

32 C. Pittman, T. Doiron, and H. Meyer, Phys. Rev. B 20, 
3678 (1979). 

33 R. D. McCarty and V. D. Arp, Adv. Cryo. Eng. 35, 1465 
(1990). 

34 B. Hetenyi, M. H. Miiser, and B. J. Berne, Phys. Rev. Lett. 
83, 4606 (1999). 




FIG. 1. Radial distribution function g(r) for 4 He at tem- 
perature T = 10 K and pressure P = 0.22746 MPa, as ob- 
tained from both classical and quantum-mechanical calcula- 
tions. The points result from simulations, while lines indicate 
second-order virial expansions. For clarity, the classical curves 
have been raised by 0.5. 
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FIG. 2. (a) Volume V per atom for 4 He (experiments cor- 
respond to the solid line, PIMD data to the circles, and the 
second-order quantum-mechanical virial expansion to the di- 
amonds) and for "classical helium" (classical MD data corre- 
spond to the open squares and the classical virial expansion 
to the dashed line) as a function of temperature at the critical 
isobar of 4 He (P cA = 0.22746 MPa). (b) The same symbols 
as in a) , but now the difference between the ideal gas volume 
Videai = ksT/P and the actual volume per atom is shown as 
a function of temperature. 
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FIG. 3. Kinetic energy (positive values) and potential en- 
ergy (negative values) for 4 He as a function of temperature 
T, at a pressure P = 0.22746 MPa. Points result from sim- 
ulations, while the solid curves result from the second-order 
virial expansion. The dashed line indicates the classical ki- 
netic energy. 



P 



25 
20 
15 
10 
5 

-5 



— P c experiment 
• P c PIMD 

- V c experiment 
o V PIMD 




16 



FIG. 4. A comparison of experimental and numerical re- 
sults for the internal energy of 4 He as a function of temper- 
ature, both along the critical isobar (solid curve and closed 
circles, respectively) and along the critical isochore (dashed 
line and open circles, respectively). In both cases, the agree- 
ment is good. 
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FIG. 5. The kinetic (positive values) and potential (nega- 
tive values) energy per atom for 4 He as a function of temper- 
ature T, along its critical isochore. Points result from sim- 
ulations (closed circles for the quantum-mechanical calcula- 
tions and open squares for the classical ones), while the solid 
curves result from second-order virial expansions. The dotted 
and the dashed curves indicate second- and third-order virial 
expansions, respectively, for the classical potential energy. 
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FIG. 6. The kinetic (positive values) and potential (neg- 
ative values) energy per atom for 3 He as a function of tem- 
perature, along its critical isochore. The closed circles were 
obtained by means of simulations and the solid curves repre- 
sent second-order virial expansions. 
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FIG. 7. Effective potential for He as obtained from the 
quantum virial expansion, for three different temperatures. 
One observes the deepening and narrowing of the potential 
well when the temperature is increased from T ~ T c to 
T ~ 4T C , 3 . The Aziz HFD-B potentiaO takes its minimum 
at r = 2.963 A. 
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